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ABSTRACT 

Gamma-ray burst afterglow flares and rebrightenings of the optical and X-ray light 
curve have been attributed to both late time inner engine activity and density changes 
in the medium surrounding the burster. To test the latter, we study the encounter 
between the relativistic blast wave from a gamma-ray burster and a stellar wind ter- 
mination shock. The blast wave is simulated using a high performance adaptive mesh 
relativistic hydrodynamics code, AMRVAC, and the synchrotron emission is analyzed 
in detail with a separate radiation code. We find no bump in the resulting light curve, 
not even for very high density jumps. Furthermore, by analyzing the contributions 
from the different shock wave regions we are able to establish that it is essential to 
resolve the blast wave structure in order to make qualitatively correct predictions on 
the observed output and that the contribution from the reverse shock region will not 
stand out, even when the magnetic field is increased in this region by repeated shocks. 
This study resolves a controversy in recent literature. 



INTRODUCTION 



Gamma-ray Burst (GRB) afterglows are produced when 
a relativistic blast wave interacts with the circumstellar 
medium around the burster and emits nontherm al radia- 
tion, (for reviews, sec Piran 2005; Meszaros 2006) The gen- 
eral shape of the resulting spectra and light curves can be 
described by combining the self-similar Blandford-McKee 
(BM) model i|Blandford fc McKed [l976T ) for a relativistic 
explosion with synchrotron radiation emission from a rela- 
tivistic electron population accelerated into a power law dis- 
tribution at the shock front. This model describes a smooth 
synchrotron light curve, with the slope of the curve a func- 
tion of the power law slope of the accelerated electrons 
and of the density structure of the surro unding medium 
jMeszaros fc Reeslll997l ; IWiiers et ai1ll997l ). 

This picture however, is far from complete and with the 
increasing quality of the available data (e.g. from swift) more 
deviations from the standard of a smoothly decaying (in the 
optical and X-ray) light curve are being found, for exam- 
ple in the shape of flares l|Burrows et aL 20051; iNousek et al.l 
l200d : lQ'Brien et alj|2006h in the X-ra y afterglows and early 
optical variability (|Stanek et al.l 120061 ) . 

Along with prolonged inner engine activity, changes 
in the surrounding density structure have often been sug- 
gested as a cause of this variability J Wang fc Loebl |2000| ; 
ILazzati et all |2002| ; INakar et ai1l200ot ). The details of the 
shape of the surrounding me dium have therefore b een the 
subject of various studies (e.g. lvan Marie et al .120061 ). as well 
as the hydrodynamics of a relativistic blast wave interacting 



with a complex density environment dMeliani fc Keppensl 
2007). Two recent studies combine a description for the 
structure of the blast wave after encountering a sudden 
change in density, like the wind termination shock of a Wolf- 
Rayet star, with an analysis of the emitte d synchrotron ra- 
diatio n that is a result of this encounter (|Nakar fc GranoH 
12001 hereafter NG and lPe'er fc Wiierj200d hereafter PW). 
but arrive at different conclusions. A short transitory feature 
in the observed light curves (at various wavelengths) is pre- 
dicted by PW, whereas NG conclude that any sudden den- 
sity change of arbitrary size will result in a smooth transi- 
tion. The purpose of this paper is to resolve this discrepancy 
in the literature by performing, for the first time, a detailed 
analysis of the radiation produced by a blast wave sim- 
ulated with a high performance adaptive-mesh refinement 
co de. For this analysis, we use the radiation code described 
in IVan Eerten fc Wiiersl (|2009l ) and th e amrvac relativistic 
magnetohydrodyna mics (RHD) code l|Keppens et al.l [20031 ; 
iMeliani et ail [2007). We take special care to perform our 
simulation at a sufficiently high spatial and temporal res- 
olution, such that a transitory feature, if any, is properly 
resolved. 

In section [2] we will first describe the setup and techni- 
cal details of our simulation run. In section [3] we will discuss 
the resulting optical light curve and the fluid profile dur- 
ing the encounter. Our numerical results confirm those of 
NG. However, by following the same approximations for the 
shock wave dynamics as PW, who approximate the different 
shocked regions by homogeneous slabs, we find that we are 
able to reproduce their result of a rebrightening of the af- 
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terglow curve. In section [4] we argue how this illustrates the 
importance of resolving the downstream density structure. 
After that we separately discuss in section [5] the contribu- 
tion of the reverse shock that is triggered when the blast 
wave hits a density discontinuity, as it is the main transi- 
tory phenomenon during the encounter. This contribution 
is overestimated by PW and assumed similar in behavior to 
that of the forward shock in NG. Since both NG and PW do 
not invoke electron cooling in their arguments and optical 
flashes, if any, occur at observer frequencies that are orders 
of magnitude below the cooling break, we will not enable 
electron cooling in our radiation code. We summarize our 
results in section [6] 



2 INITIAL CONDITIONS 

We will study the case of a massive (M > 25M©), low metal- 
licity (Z ^* 0.01/70 ) progenitor star. During its Wolf-Rayet 
phase (lasting ^ 10 6 years) a stellar wind is produced, which 
determines the shape of the circumstellar medium. The typ- 
ical mass-loss rate is approximately M ^ 10 _6 Mq yr _1 
and the typical wind velocity v w ^ 1000 km s _1 . Because 
the stellar wind flow is supersonic, a shock is produced. A 
simple schematic description of the circumstellar medium 
(where we ignore complications such as the influence of 
photo-ionization) consists therefore of (starting near the star 
and moving outwards) a free-flowing stellar wind region, a 
density jump separating the stellar wind region from a ho- 
mogenized region influenced by the reverse shock, a contact 
discontinuity followed by a region shocked by the forward 
shock. The forward shock front then separates the shocked 
medium from the unshocked interstellar medium (ISM). 

Following the GRB explosion, a relativistic blast wave is 
sent into this environment. For the typical progenitor values 
above, an ISM number density uism ^ 10 3 cm -3 and a 
GRB explosion energy of E = 10 53 erg, this blast wave will 
only encounter the first discontinuity during its relativistic 
stage. The discontinuity will be positioned at Ro — 1.6 • 
10 18 cm and corresponds to a jump in density of a factor 
4. Before the jump the radial density profile is given by 
n(r) — 3 • (r/1 ■ 10 17 )~ 2 cm" 1 , and after the jump by the 
constant n(r) =4x3- (Ro/1 ■ 10 17 ) -2 cm" 1 . These exact 
values are chosen to conform to PW. 

We have run a number of simulations of relativistic blast 
waves hitting the wind termination shock at Rq. The initial 
fluid profile is generated from the impulsive energy injection 
BM solution with the parameters described above for the 
explosion energy and circumburst density, keeping the adi- 
abatic index fixed at 4/3. The starting time is taken when 
the shock Lorentz factor is 23. The blast wave will hit the 
discontinuity when its Lorentz factor is ^ 22.27, at an ex- 
plosion lab frame time t enc = 5.34 • 10 7 seconds (with t = 
set to the start of the explosion). This time corresponds to 
^ 0.3 days for radiation coming from the shock front in ob- 
server time (which is taken to be zero at the start of the 
explosion). To completely simulate the encounter, we will 
follow the evolution of the blast wave from 5 ■ 10 6 seconds 
to 6.4 • 10 7 seconds and will store enough output to obtain 
a temporal resolution (in lab frame simulation time) dt of 
1.56 • 10 3 seconds. 

For the outer boundary of the computational grid we 



take 6-10 cm, enough to completely capture the shock 
profile during the encounter even if it were to continue at 
the speed of light. In order to resolve the shock wave, even 
at its smallest width at Lorentz factor 23, we take 10 base 
level cells and allow the adaptive mesh refinement routine to 
locally double the resolution (where needed) up to 17 times. 
This implies an effective resolution dr ^ 6.3 ■ 10 11 cm and 
effectively 1,310,720 grid cells. 

Three simulations were performed using the initial con- 
ditions from PW (along with some at lower resolutions, to 
check for convergence): a test run with stellar wind pro- 
file only (and no discontinuity), one with a density jump 
of 4 and one with a far stronger density jump of 100. Al- 
thoug h density jumps muc h larger than 4 may be feasible 
(see Ivan Marie et alj 120061 . for an example scenario where 
the progenitor star has a strong proper motion -the relativis- 
tic blast wave will then be emitted into a stellar environment 
that takes the shape of a bow shock), this is not the main 
motivation for the factor 100 simulation run. The primary 
focus is on establishing if the lack of an observer effect in the 
light curve persists for general values of the density jump. 

To study relativistic as well as ultra-relativistic blast 
waves, in addition to the Lorentz factor 23 scenario we have 
also performed two simulations (one with jump and a test 
run without) where we moved the density jump outward to 
3 ■ 10 19 cm, while keeping the other parameters equal. In this 
scenario the blast wave encounters the jump when it has a 
shock Lorentz factor ^ 5. 

The simulation output is then analyzed using the ra- 
diation code for an observer at a distance of 1 • 10 28 cm. 
The microphysics of the shock acceleration is captured by 
a number of ignorance parameters. The fraction of thermal 
energy residing in the small scale downstream magnetic field 
is e_B = 0.01, the fraction of thermal energy in the acceler- 
ated particle distribution e_B = 0.1, the number of power 
law accelerated electrons as a fraction of the electron num- 
ber density £jv = 1 and the slope of this power law p — 2.5. 
Again these values are chosen to match PW. 



3 LIGHT CURVE AND SHOCK PROFILE 

The discussion below refers to the shock Lorentz factor 23 
scenario. The Lorentz factor 5 simulations lead to qualita- 
tively similar light curves and will therefore not be discussed 
in further detail. The transition then takes extremely long 
due to the longer dominance of earlier emission. These simu- 
lations confirm that the results presented hold for relativistic 
blast waves as well, not just for ultra-relativistic blastwaves. 

Directly after hitting the discontinuity, the blast wave 
splits into three regions. The innermost region, up to the 
reverse shock (RS) front remains unaware of the collision. 
Beyond the RS the plasma gets homogenized up to the con- 
tact discontinuity (CD). The region following the contact 
discontinuity, up to the forward shock (FS) is not homoge- 
neous but will gradually evolve into a BM profile again for a 
modified value of the circumburst density structure. A snap- 
shot of the shock structure during the encounter is shown in 
figure [T] We show comoving density (as opposed to the lab 
frame density) because the differences between the different 
regions then stand out more clearly. 

The optical light curves calculated from the simulations 
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Figure 1. A snapshot of the comoving density profile at 17 re- 
finement levels of the fluid at emission time t e = 5.48578 • 10 7 
s, for the factor 4 increase in density. The different regions are 
clearly visible. From left to right we have: up to the steep rise the 
region not yet influenced by the encounter, the plateau resulting 
from the passage of the reverse shock, and starting at the gradual 
rise the region of the forward shock. The front part of the for- 
ward shock region is again homogeneous in density, showing the 
difference between the idealized BM solution and actual simula- 
tion results. The flat part of the forward shock region (smallest, 
rightmost region) is resolved by ^ 100 cells. 



are observed at v = 5 ■ 10 14 Hz, which lies between the syn- 
chrotron peak frequency v m and the cooling break frequency 
v c (it may be helpful to emphasize that here, contrary to 
shock interaction during the prompt emission phase, v m is 
found at a similar frequency for both the forward and re- 
verse shock contributions). Because the observer frequency 
lies well below the cooling break, we ignore the effect of elec- 
tron cooling. The light curves for the factor 4 and factor 100 
density jumps are found in figure[2] For complete coverage at 
the observed times and clarity of presentation, analytically 
calculated emission from a BM profile with Lorentz factors 
> 23 (or > 5) has been added to that calculated from the 
simulations. From the light curves we draw the following 
conclusion: an encounter between the relativistic blast wave 
and a wind termination shock does not lead to a bump in 
the light curve, but instead to a smooth change in slope. The 
new slope eventually matches that of a BM solution for the 
density structure found beyond the discontinuity. 



4 RESOLVED BLAST WAVE VERSUS 
HOMOGENEOUS SLAB 

The optical light curves presented in the above section differ 
distinctly from those presented by PW in that they show no 
bumps. This difference in results has to be caused by one or 
more differences in our assumptions, which are: 

• PW include both electron cooling and synchrotron self- 
absorption, while in this paper we have included neither. 

• We take the magnetic field to be a fixed fraction es of 
the local thermal energy in all parts of the fluid, even those 
shocked twice, whereas PW have a magnetic field in the 
reverse shock region that is slightly higher. This is because 
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Figure 2. The figure shows the resulting optical light curves at 
5-10 14 Hz, for the cases of a continuous stellar wind environment, 
a jump of a factor 4 followed by a homogeneous environment and 
a jump of a factor 100. 50 data points have been devoted to 0.3 
- 1 day and 50 data points to the following 19 days. A smooth 
transition towards the power law behaviour corresponding to a 
BM shock wave expanding into a homogeneous environment is 
visible, even for the extreme change in density. 



analytical approximation of dynamics 
same approximation w/o advected B-1ield 
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Figure 3. Resulting light curves at 5- 10 14 Hz when our radiation 
code is applied to the homogeneous slabs approximation of PW, 
instead of a hydrodynamical simulation. The bottom curve shows 
the resulting light curve if the magnetic field in the reverse shock 
region does not contain the additional increase in magnetic field 
in the reverse shock region. Contrary to the light curves shown in 
figure [2] in both cases a clear rise in intensity with respect to the 
previous level is seen over the course of a few hours, as predicted 
by PW for homogeneous slabs. 



they take into account that the dominant magnetic field in 
the reverse shock region is actually the field advected with 
the flow from the region shocked once. The newly created 
field is approximately a factor 1.2 smaller. 

• We resolve the downstream fluid profile, while PW ap- 
proximate the different regions behind the shock front by 
homogeneous slabs of varying density, thermal energy and 
Lorentz factor. Also, they freeze the fluid Lorentz factors 
during the encounter. 
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Since the optical light curve corresponds to an observer 
frequency sufficiently above the self-absorption critical fre- 
quency and sufficiently below the cooling break frequency, 
neither cooling nor absorption should have any visible ef- 
fect on the shape of the curve. The fact that cooling is not 
required for the bump found by PW is also immediately 
obvious from figure [3] where we have applied our radiation 
code directly to the homogeneous slabs approximation of 
PW, with electron cooling disabled. The light curves thus 
generated do show a bump feature after the onset of the 
encounter (this also provides a check on the internal consis- 
tency of both models). To explicitly check the effect of the 
stronger magnetic field in the reverse shock region we have 
generated two light curves: one where all the fluid quantities 
are exactly similar to those of PW and one where we ignored 
the stronger field in the reverse shock region but kept the 
field at fixed fraction of the thermal energy (which is the 
same as that in the forward region in the homogeneous slab 
approximations, due to pressure balance across the contact 
discontinuity). As can be seen from the figure, the tempo- 
rary rise occurs in both cases, with only a marginal difference 
between the two curves. 

This brings us to the third difference listed. We conclude 
that to determine the visible response of a blast wave to den- 
sity perturbations, it is crucial to take the radial structure of 
the blast wave into account. This (along with establishing 
the lack of a transitory feature itself) is the main conclusion 
from this paper and forms an important justification for the 
kind of detailed approach that we have employed, where the 
dynamics of the blast wave are simulated using a high per- 
formance RHD code, together with a radiation code that 
accurately probes all local contributions to the synchrotron 
spectrum. It is also important to emphasize that the bump 
found by PW is not the result of inaccurately modeling the 
different arrival times for photons arriving from different an- 
gles relative to the line of sight, as has been stated in NG. 
This also can be seen from figure [3] which confirms that, for 
homogeneous slabs, the light curves published by PW are 
calculated correctly. 

The importance of the downstream shock structure can 
be understood as follows. By taking a homogeneous slab one 
not only locally overestimates the downstream density, but 
also the Lorentz factor and thermal energy (and hence the 
magnetic field). Also, the width of the homogeneous slab 
is determined by comparison with the downstream density 
structure or the energy density structure or the velocity 
structure, and matching the width to one of these comes at 
the expense of a lack of similarity to the others. (And finally, 
keeping the Lorentz factors fixed during the encounter also 
contributes to the overestimation of the flux emitted during 
the encounter). Essentially, all this indicates a lack of reso- 
lution. The homogeneous slab implies a spatial resolutiorQ 

1 Even though PW identify three different regions during the 
encounter, this in itself does not imply an improved spatial res- 
olution, since the fluid conditions in each region are connected 
to each other (and the upstream medium) via shock-jump con- 
ditions that strictly speaking require all regions to be directly 
adjacent at the same position. The simulation snapshot in fig. 
[T] shows that the assumption of the reverse shock region being 
thermalized and isotropic is not unreasonable, but also shows a 
clear density gradient within the forward shock region. 
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Figure 4. Received flux at observer frequency 5 ■ 10 Hz, cal- 
culated for a single emission time t e = 5.48578 ■ 10 7 s (the same 
time as in fig. ffj. Curves are shown both for the homogeneous 
slab approximation and for the numerical simulation fluid pro- 
file. In each case the contribution from the different regions has 
been marked: the top curve shows the total flux, the curve below 
the flux when the contribution from the forward shock region is 
omitted and for the lowest curve the reverse shock region has been 
omitted as well. The flux level for the homogeneous slab approx- 
imation is much higher than that from the simulation, with (for 
this particular emission time) the contribution from the reverse 
shock dominating the total output. At the same emission time, 
the reverse shock region contribution for the simulation is still 
significant, but no longer dominant. For the simulation snapshot 
we have estimated the position of the contact discontinuity, and 
therefore the edge of the reverse shock region, at the right edge 
of the plateau, before the onset of the rise in density, (see fig.[T]) 

Ar ^ R/T 2 cm (with R the blast wave radius and F the 
blast wave Lorentz factor) , and is therefore in principle only 
applicable to describe behavior on time scales At > Ar/c. 
This is true in general, not just for simulations, and in our 
case yields At ^ 1.5 days at the time of the onset of the en- 
counter. The reason that the homogeneous slab does work 
to describe the general shap e of the light cur ve from the BM 
blast wave, as was done by IWaxmanl (|l997f ) among others, 
is that in these cases the slab is used to describe behavior 
on time scales At >> Ar/c (actually At arbitrary large, 
for understanding of asymptotic behavior). But one should 
for example not expect the homogeneous slab approxima- 
tion to get the absolute scale right, and indeed it is off by 
a factor of a few (justifying more detailed calcu lations like 
iGranot fc Sar1l2002l : IVan Eerten fc Wiiedl2009h . 

5 THE REVERSE SHOCK CONTRIBUTION 

In the previous we have established that the reverse shock 
caused by the encounter with the density perturbation does 
not cause a rise in the observed light curve. Since this re- 
ve rse shock has been evoked to explain rebrightening (e.g. 
by iRamirez-Ruiz et al.l [20051 ). it is of interest to look at its 
contribution in some more detail. In fig. [4] this contribution 
(in the optical) is compared directly to the total flux emit- 
ted from the shock profile, both for the simulation and for 
the PW approximation. The important difference is the rel- 
ative overestimation of the reverse shock region in the PW 
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approximation. The relative contributions for the different 
regions within either the homogeneous slab or the resolved 
blast wave simulation of course depend on their relative sizes 
and therefore on the emission time. Another feature of note 
is that the homogeneous slabs approximation results in an 
emission profile that is sharply peaked, whereas the more 
accurate profile displays a flatter tail and a smoother tran- 
sition between rise and decay. 

The shock structure is also calculated and implemented 
in NG, starting from the shock jump conditions and assum- 
ing homogeneous slabs for the forward and reverse shock re- 
gion, yet they do not find a temporary rebrightening. This 
is a consequence of the fact that they set the reverse shock 
contribution at a fixed fraction of the forward shock con- 
tribution, while allowing this forward shock contribution to 
evolve according to the appropriate BM profile following the 
density change, as opposed to freezing the shock Lorentz 
factors during the encounter. That the forward shock de- 
termines the shape of the light curve is then imposed as a 
feature of their model (i.e. in their equation 20) and yields 
an adequate heuristic description of the light curve found as 
a result of their simulations. 

The difference between the simulations by NG and ours 
is merely a technical one: instead of a Eulerian code (that 
can also be used for simulations in more than one dimension, 
which we will perform in future work), they use a Lagrangian 
code for the dynamics. The reconstruction of the light curves 
from the code is equivalent. They also, like us, do not take a 
slight increase in the magnetic field in the reverse shock re- 
gion into account. NG provide no information on the spatial 
and temporal resolution of their simulations. 



6 SUMMARY AND CONCLUSIONS 

We have performed high resolution hydrodynamical simula- 
tions of a relativistic blast wave encountering a wind termi- 
nation shock and have calculated th e resulting light curve 
using the radiation code described in lVan Eerten fc Wiiersl 
(2009). As a result we have found no variability in the 
optical, not even for very large density changes, for blast 
waves in the self-similar phase. This renders it very un- 
likely that observed optical variability in GRB afterglow 
light curves can be explained from density perturbations 
in the externa l medium surroundin g the burster, as sug - 
gested by e g IWang fc Loebl (|200Ch : lLazzati et ail ||2002n ; 
iNakar et al l (|2003h . PW. This research, however, has been 
limited to spherically symmetric density perturbations. A 
second caveat is the assumption of self-similarity for the 
blast wave appr oaching the wind t e rmina tion shock. As 
demonstrated bv lMeliani fc Keppensl (|2007) , for a termina- 



tion shock close to the star (R ^ 10 cm in their simulation, 
for a short Wolf-Rayet phase), the blast wave structure may 
still somewhat retain the initial structure of the ejecta (in 
their simulation, a uniform static and hot shell, i.e. fireball), 
which may have observable consequences. The latter is how- 
ever not likely, given the already reasonably strong resem- 
blance between their simulation output during the encounter 
and ours, where the same shock regions can be identified in 
the fluid profile with similar values for the physical quan- 
tities of interest. Also, if the pre-encounter shock wave is 
sufficiently different from the self-similar solution this will 



also have consequences for the global shape and temporal 
evolution of the observable light curve, and the slope will 
become markedly different from the one predicted from the 
BM solution. 

Of the two main explanations for (sometimes quite 
strong) late optical variability, refreshed or multiple shocks 
appear to be a far more realistic option than circumburst 
medium interactions. We are currently performing simula- 
tions on multiple interacting shocks to test this alternative 
hypothesis. 

We have compared the results of our simulation to the 
literature and from a comparison to the approximations and 
assumptions used by PW and NG especially, we conclude 
that the fact that we resolve the radial blast wave structure 
explains the discrepancy between our resuls and those of 
PW. This, in turn, forms an important justification for the 
kind of detailed approach that we have employed, where the 
dynamics of the blast wave are simulated using a high perfo- 
mance RHD code, together with a radiation code that accu- 
rately probes all local contributions to the synchrotron spec- 
trum. We note that, contrary to what is stated by NG, the 
calculation of angula r smearing of th e signal in PW (which 



in turn was based on I Waxmanll f 9971 ) is correct 
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